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CHAPTER  I 


INTRODUCTION 

In  applications  of  the  theory  of  random  processes,  a function 
which  is  often  important  is  the  probability  density  function  of  the 
duration  t,  of  surges  of  a random  process  above  a given  critical 
value.  A practical  method  for  deriving  this  probability  density 
p^(a),  for  a given  random  process  remains  \inknown.  Thus  the  engineer 
who  needs  to  know  p^(a)  must  be  content,  instead,  with  approximations 
to  it.^  The  approximation  most  widely  accepted  was  published  in 
1958  by  S.  0.  Rice.^  This  approximation  has  been  proved  to  be 
accurate,  but  it  is  difficult  to  obtain.  In  addition.  Rice's  work 
only  applies  directly  to  Gaussian  processes  and  is  different  to 
extend  to  non-Gaussian  processes. 

A.  N.  Denisenko  has  recently  proposed  a new  approximation  which 
is  far  easier  to  evaluate  for  Gaussian  processes  and  which  can  be 


^The  notation  p(t)  is  often  used  for  this  distribution,  but  to 
avoid  confusion  between  the  random  variable  and  the  argument  of  its 
distribution  function,  p^(a)  will  be  used  in  this  paper.  With  this 
notation,  the  probability  of  t being  less  than  x is 

f P^(a)do( 

J —00 

^Stephan  0.  Rice.  "Distribution  of  the  Duration  of  Fades  in 
Radio  Transmission:  Gaussian  Noise  Model."  The  Bell  System 
Technical  Journal.  Vol.  37,  No.  3 (May  1958),  pp.  581-635. 
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easily  extended  to  non-Gaussian  processes.^  In  his  presentation 
of  this  new  approximation,  Denisenko  included  computer  simulation 
data  which  tended  to  verify  the  accuracy  of  his  approximation. 

The  p\irpose  of  this  thesis  is  to  check  the  results  of 
Denisenko's  approximation  against  experimental  data,  including 
non-Gaussian  processes. 


^A.  N.  Denisenko.  "Estimate  of  the  Distribution  of  Surge 
Durations  for  Random  Processes."  Radiotekhn , i Elektr . Vol.  20, 
(July  1975),  PP.  1529-1532. 


CHAPTER  II 


THEORETICAL  BACKGROUND 
The  Mean  of  the  Surge  Distribution 


While  the  exact  solution  for  p^(a)  remains  unknown,  the  mean  of 
the  distribution  is  known.  For  a random  process  X(t),  which  at  any 
given  time  has  probability  distribution  “the  mean  duration  of 

a surge  above  level  x^  is 


E[x] 


P(X>x  ) 
o 


where 


P(X>x  ) 
o 


r p^(e)  dg 

o 


(1) 


and  V is  the  average  frequency  with  which  X crosses  x^  given  by 


v{xj  = |^JP^^^,(x^,y)  dY 


where  p^  ^,(a,b)  is  the  two  dimensional  probability  distribution  of 
X(t)  and  X’(t).l 


^Petr  Beckmann.  Probability  in  Communication  Engineering. 

New  York , Harccrt , Brace~&~WorldIncT7~l9^7T'p'^~2297*'Err''isthe 
expected  value  operator,  and  a prime  indicates  the  derivative  with 
respect  to  the  argiiment. 


For  a normal  process  having  mean  y,  variance  a , correlation 
function  B(t)  = E[X(t)  X(t+T)],  and  correlation  coefficient  R(t) 

= B(t)/o  ; equation  (l)  gives^ 


E[t]  = 


(x  -y)  X -y 

o _ o 

exp  — erfc 


/-R"(0)  2a 


where 


erfc(z)  = 1 - erf(z)  = 1 - — 


z 2 
-o  , 
e da 


(2) 


If  the  process  X(t)  is  stationary,  exponential,  "analogus  to  a 
normal  process",  and  with  derivative  X'{t)  independent  of  X(t),  then 
R(t)  = v^B(t)  and^ 


(3) 


It  is  Interesting  to  note  that  this  mean  value  is  independent  of  the 
surge  level  x^- 

While  these  results  are  important  and  useful,  p^(a)  must  still 
be  estimated  because  the  above  calculations  yield  only  the  mean  of 
the  distribution  and  say  nothing  about  how  p^(ci)  is  distributed  about 
the  mean . 


2lbid.  p.  230, 

^Petr  Beckmann,  Orthogonal  Polynomials  for  Engineers  and 
Physicists,  Boulder , Colorado,  The  Golem  Press,  1973,  p,  196, 

The  expression,  "analogus  to  a normal  process",  means  that  the  two 
dimensional  density  of  the  process  can  be  expanded  in  orthogonal 
polynomials  in  a cononical  form  as  explained  in  this  reference. 


Rice's  Approximation 

The  approximation  to  the  surge  distribution  which  Rice  pro- 
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‘'Stephan  0.  Rice.  "Distribution  of  the  Duration  of  Fades  in 
Radio  Transmission:  Gaussian  Noise  Model."  The  Bell  System 
Technical  Journal,  Vol.  XXXVIII,  No.  3 (May  1958)»  P-  _601. 
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2(l-r^) 


Next  Rice  approximates  the  distribution  of  very  long  surges 
(laj-ge  values  of  a where  pj^^(a)  is  no  longer  even  closed  to  valid) 


These  two  approximations  (equations  (U)  and  (5))  are  then 
plotted  and  a value  of  a is  picked  below  which  (4)  will  apply  and 
above  which  (5)  will  hold.  At  the  same  time  and  are  adjusted 
so  that  ( i ) the  resulting  curve  looks  "reasonable" , ( ii ) the 
resulting  approximation  p^(a)  has  the  proper  mean  (given  by  equation 


(2)),  and  finally  ( 


iii)  p , 


)da  = 1. 


All  of  this  is  done  by  eye. 


so  we  see  that  artwork  and  experience  as  well  as  mathematics  play  a 
role  in  obtaining  p^(a)  by  Rice's  method.  It  has  been  amply  verified 
that  approximations  obtained  by  Rice's  method  agree  well  with  experi- 
mental data.^ 


Denisenko's  Approximation 

Denisenko  approached  the  problem  by  estimating  the  cumulative 


distribution  function  F ( 

T 


^ r\ 


(.y)dy 


which  is  equal  to  the 


Sibid.  pp.  589-599. 


0 


probability  of  a sxirge  of  length  t^o,  which  is  equal  to  1 minus  the 
probability  of  a surge  of  length  T>a.  The  exact  expression  for  the 
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probability  of  a surge  of  length  t>o  can  be  expressed  as® 

lim[P(X(t_)>x  , X(t  )>x  ,...X(t  )>x  |x(t  )>x  )]  (6) 

n4co  2 o 3 o n o 1 o 

where  t -t,  = a and  the  t are  all  chosen  so  that  t , -t  = a/n. 

H X K ic  ■ X A 

Let  U,  be  the  event  X(t,  )>x  . Using  this  to  rewrite  equation  (6) 
k k o 

gives  the  cumulative  distribution 

F^(a)  = 1 - lm[P(U2|u^)'P(U2|U^,U2)...P(U^|U^.U2,U^....U^_^)] 

(7) 

At  this  point , Denisenko  states  without  proof  that 

P(UJUi,U2,U3....U^.^)<P(U^lU^_l)  (8) 

He  then  uses  the  lower  bound  of  this  expression  as  an  approximation 
to  the  left  side,  ie. 

P(UJUi.U,,U3....U^.l)  ■P(UjU^.i)  (9) 

Equation  (9)  applied  n-2  times  in  equation  (7)  gives 

P,(a)  = 1 - lialPtUnlVl”” 


or 


F^(a)  = 1 - lm[P(X(tj^)  > x^lX(tj^_^)  > x^)]" 


®P(  ) is  the  probability  of  the  argument  in  brackets,  a vertical 
bar  means  "given  that",  and  a comma  means  "and". 


F^(o)  = 1 - l^[P(X(t+^)  > x^lx(t)  > x^)]" 


and  finally,  Bayes  Rule  yields 

^ P(X(t+— ) > X , X(t)  > X ) 

F (a)  = 1 - 1^[ ^ 2 2_]  (10) 

^ P(X(t)  > X ) 

o 

Up  to  this  point,  the  density  P^(o)  of  the  process  X(t)  has  not 

been  considered.  Now  equation  (lO),  when  applied  to  a normal 

2 

process  with  mean  \i  and  variance  a gives  (following  a Taylor  series 
expsinsion  about  R(0)  and  some  algebra) 


F (a)  = 1-e 

T 


-Aa 


where 


2Tr[l-<t— 2 ] 


2a 


and 


Hz)  = -= 


g-t^/2  dt 


(11) 


(12) 


(13) 


And  finally  taking  the  derivative  of  equation  (ll)  with  respect  to 
o,  yields  the  estimate  of  the  density  of  t. 

p^(a)  = Ae”'^“  (lU) 


For  an  unbiased  approximation,  the  mean  of  this  estimate 
should  agree  with  the  known  mean  of  the  surge  length  distribution. 
As  shown  in  Appendix  E,  this  is  the  case. 


V ^ 


The  Extension  of  Denisenko's  Approximation 
The  direct  application  of  equation  (lO)  to  non-normal 
processes  results  in  expressions  which  cannot  be  evaluated  in 
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closed  form.  However,  working  directly  from  equation  (12)  Beckmann 
extended  Denisenko's  resiilts  to  non-normal  processes.^ 

Suppose  that  the  surge  distribution  of  a process  Y(t)  is 
required,  where  Y(t)  has  arbitrary  probability  distribution  Py(a) 
and  corresponding  cumulative  distribution  function  FyCo) 


j>^{y)dy. 


Let  the  process  X(t)  be  normally  distributed  with 


zero  mean  and  unit  variance.  Thus  X(t)  has  cumulative  distribution 
function  $(o)  given  by  equation  (13). 

The  process  Y(t)  is  obtained  as  some  monotonic  function  G(X) 
of  X(t),  and  G(X)  is  picked  so  that  Y(t)  has  the  desired  distribution 
function. 


Y(t)  = G(X(t)) 

P(X  £ a)  = P{Y  j<  G(a)) 

$(a)  = Fy(g(o)) 

G(a)  = Fy'^CtCa))  (15) 


^Petr  Beckmann.  "Probability  Distribution  of  Surges  and 
Fades."  Proceedings  of  the  IEEE,  Vol.  6I4,  No.  U,  (April  1976), 
pp.  571-572. 
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In  this  way  it  is  always  possible  to  find  a function  G which 
transforms  a normal  distribution  to  the  desired  distribution  of  Y. 
The  duration  of  a surge  of  X(t)  above  level  will  be  equal  to  the 
duration  of  a surge  of  y(t)  above  level  G(x^).  Thus  to  get  the 
surge  distribution  of  Y(t)  above  level  y^,  equation  (ll)  can  be 
used  with  the  terms  in  equation  (l2)  replaced  by  G 
is  important  to  remember  that  the  term  /-R"(0)  in  equation  (12) 
now  refers  to  the  process  Y(t)  and  Ry(t)  must  be  evaluated  from 


R (r) 
y 


fCO 


—00 


r G(x.)G(x_)^ 
[-  1-  -2-3  . 

• 2ir A-R^(t) 

X 


X -x^ 

exp[-^ ^ 1 g-ldx.  dx. 

2(1-R^(t))  1 2 


(16) 


where  a is  the  variance  of  Y(t). 

y 

The  evaluation  of  equation  (l6)  can  be  simplified  somewhat  by 
use  of  orthogonal  polynomials.  However,  even  if  orthogonal  poly- 
nomials are  used,  evaluating  this  equation  (since  the  second 
derivative  of  the  result  is  required),  and  finding  G(a)  can  involve 
difficult  computations.  Fortunately  a simpler  method  exists  which 
can  eliminate  these  ciombersome  calctilations . 

The  mean  of  the  surge  distribution  of  Y(t)  can  be  calculated 
directly  from  a knowledge  of  Y(t)  using  equation  (l).  Because  the 
mean  of  Denisenko's  approximation,  and  the  mean  of  the  true 
surge  distribution  coincide  exactly,  the  reciprocal  of  the  mean 
obtained  from  equation  (l)  can  be  substituted  for  A in  equation  (l4) 


and  the  problem  is  solved. 


For  example,  if  the  surge  distribution  of  an  exponential 
process  is  desired,  examination  of  equations  (3)  and  (lU)  effort- 
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lessly  give 


p^(a)  = 


-R”(0) 

8 


exp[- 


-R"(0) 

8 


a] 


(17) 


It  is  interesting  to  note  that  this  approximation  to  the  surge 

distribution  of  an  exponential  process  is  independent  of  the  surge 

level  X . 
o 

Characteristics  of  the  New  Approximation 
The  key  to  understanding  Denisenko's  approximation  to  the  surge 
distribution  lies  in  understanding  the  implications  of  equations 
(8)  and  (9).  as  these  are  the  only  approximations  made  in  the 
derivation  of  the  estimate,  equation  (lU),  from  the  exact  expression, 
eqiiation  (6). 

In  Denisenko's  paper  one  of  the  basic  premises,  rewritten  here 
as  equation  (8),  was  given  without  proof. 

GDI  _>  CD2 

where 

CDl  = conditional  density  1 = P(X(t,  ) > x |x(t  ) > x ) 

K O it“X  O 

and 

CD2  = conditional  density  2 

= P(X(t^)  > X |x(t, ) > X ,X(t^)  > X ,...X(t.  ,)  > X ) 
k o 1 o 2 o k— 1 o 

Critical  scrutiny  shows  this  claim  to  be  false.  A counter- 
example is  outlined  in  Appendix  A where  it  is  shown  that  CDl  can 
also  be  less  than  CD2.  This  leaves  no  particular  magnitude 


J 
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I 

1 relationship  "between  CDl  and  CD2.  If  equation  (8)  were  true,  then 

^ X ^ « 

the  approximation  of  equation  (9)  CDl  = CD2  would  introduce  a 
consistent  bias  into  the  estimate  of  the  surge  distribution,  i.e. 
the  estimated  cumulative  distribution  function  would  always  be 
greater  than  the  actual  cumulative  distribution  function.  Because  i 

equation  (8)  fails,  no  consistent  bias  should  be  expected  in  the 

I estimate  of  the  stirge  distribution. 

: !■ 

! 

P(x(o.  ) > x^lx(t-)  > x^,X(t_)  > x^,...X(t,  ) > X ) 

K O 1 O ti  O iC-1  O 

If  all  of  the  greater  than  signs  in  the  conditions  of  equation 
(9)  are  replaced  with  equal  signs,  then  equation  (9)  would  be  true 
for  a Markov  process.  If  X(t)  is  a normal  process  and  Markov,  then 
by  Doob's  Theorem,  X(t)  has  an  exponential  correlation  function, 

R(t).®  Thus  if  Denisenko's  estimate  is  used  to  approximate  the 
stirge  distribution  of  a normal  process  with  exponential  correlation 
function,  an  estimate  which  is  reasonably  close  to  the  true  distri- 
bution would  be  expected  and  inaccuracies  in  the  approximation  would 
be  due  to  the  greater  than  signs  instead  of  equal  signs  in  equation 
(9)-  If  Beckmann's  extension  is  used  for  non-normal  processes  with 
exiMDnential  correlation  fiinctions,  further  inaccuracies  should  be 

®J.  L.  Doob.  Stochastic  Processes.  New  York,  John  Wiley  & 

Sons,  Inc.,  1953,  p.  233.  This  result,  commonly  known  as  "Doob's 
Theorem",  was  never  stated  as  a theorem  by  Doob.  He  gave  this 
result  for  normal  processes  as  a special  case  of  a much  more  general 
theorem  on  correlation  functions  of  Markov  processes. 


k 
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expected  "because  Doob's  Theorem  no  longer  applies  for  the  non-normal 
cases,  unless  they  are  analogus  to  a normal  process  (see  footnote  3 
of  this  chapter).  If  a process  Y(t)  is  obtained  as  a non-linear 
mapping  of  a normal  process,  as  was  done  in  this  case,  then  y(t) 
is  not  analogus  to  normal. 

* ! 

I i If  the  approximation  is  applied  to  a process  with  arbitrary 

I ! correlation  function,  the  accioracy  of  the  approximation  would  be 

j very  sensitive  to  how  well  equation  (9)  holds  for  the  particular 

correlation  function  in  question. 


CHAPTER  III 


EXPERIMENTAL  PROCEDURE 

To  exajnine  Denisenko's  approximation  to  the  surge  distrihution, 
a computer  program  was  used  which  models  the  process  X(t)  by  a 
digital  Monte  Carlo  simulation  and  records  the  duration  of  surges 
above  the  level  x^.  The  program  is  capable  of  simulating  X{t)  as 
either  a normal  process  or  an  exponential  process  which  allows 
evaluation  of  Beckmann's  extension  to  non-Gaussian  processes  as  well 
as  the  estimate  of  Gaussian  processes. 

The  Second  Order  Autoregressive  Process 
As  indicated  at  the  end  of  Chapter  II,  it  is  necessary  to 
evaluate  the  performance  of  the  estimators  for  processes  with  many 
different  types  of  correlation  functions.  To  give  the  program  this 
flexibility,  X(t)  was  generated  using  a second  order  autoregressive 
process.  The  discrete  second  order  autoregressive  process  with  zero 
mean  is 

where  is  independent  of  X^^  for  all  k,  and  is  white.  Adjust- 

ments of  the  coefficients  A1  and  A2  in  equation  (l8)  make  possible  a 
wide  reinge  of  autocorrelation  functions  of  Xj^. 
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Generating  the  Gaussian  Process 
Because  the  second  order  autoregressive  process  is  linear,  an 
input  process  Z which  is  Gaussian  will  result  in  an  output  process 

iv 


which  is  also  Gauss i£in 

"direct  method"  was  used^ 


To  achieve  the  uncorrelated  normal  process  the  so-called 


Z,  = /-P  In  U,  cos  (PttU,  ) 
k k k+1 


where  the  varieties  are  independent,  uncorrelated  and  uniformly 
distributed  on  the  interval  from  zero  to  one.  The  values  of  Uj^  were 
obtained  from  the  function  RAUF  of  the  CDC  FORTRAN  library.  This 
gives  Z with  mean  zero  and  variance  one.  Because  the  output  of  the 
autoregressive  process  is  also  required  to  have  zero  mean  and  unit 
variance,  the  output  must  be  scaled  as  shown  in  Appendix  B. 


Beginning  in  the  Steady  State 

The  process  generator  which  results  from  the  considerations  of 
the  previous  two  sections  is 


with  K the  normalization  constant  given  in  Appendix  B. 


^Milton  Abramowitz  and  Irene  Stegun,  Handbook  of  Mathematical 
Functions , National  Bureau  of  Standards.,  196U,  p.  953. 
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V/hen  the  simulation  is  started,  initial  values  for  ^ and 
2 must  be  chosen  so  that  the  process  is  in  the  steady  state  from 
the  first  sample  onward.^  When  the  generator  is  in  the  steady 
state,  the  values  Xj^  ^ and  X^^  ^ will  be  normally  distributed  with 
zero  mean,  unit  variance,  and  a non-zero  covariance,  R. 

In  Appendix  C it  is  shown  that  the  proper  choice  of  X*  ^ and 
X*  g,  where  the  astrisks  indicate  initial  values,  is  given  by 


^-1  = "l 

and 

X«_2  = KZi  * 


where  the  values  and  Zg  are  independent  samples  drawn  from  a 
distribution  which  is  normal  with  zero  mean  and  unit  variance. 


Generating  the  Exponential  Process 
The  exponential  process  was  generated  by  taking  a fxinction  G 
of  the  normal  process  outlined  in  the  previous  sections. 

Let  the  process  Y be  exponentially  distributed 

Py(y)  = J exp[-  for  o _<  y < “ 

= 0 for  y < OL 


^The  process  generator  being  in  the  steady  state  means  that 
the  process  is  stationary. 
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30  that  the  mean  of  Y is  a+3  and  the  variance  of  Y is  3 . The 
cumulative  distribution  of  Y is 

P(Y  £ y)  = 1 - exp[-  for  a ^ Y * 

= 0 for  y < a 

The  cumulative  distribution  function  of  the  normal  process  X is  the 
♦ function  of  equation  (13).  The  function  G,  computed  from  (l5),  is 


G{x)  = 3[-  ln(^  erfc  (l^))  + a] 


(19) 


This  solves  the  problem  of  mapping  the  normal  process  X into  an 
exponential  process  Y.  The  computation  of  this  function  is 
accomplished  as  shown  in  Appendix  D. 


The  Process  Probability  Density  Function 
To  give  a quick  visual  check  of  the  density  of  the  process 
X(t),  a histogram  of  X(t)  is  plotted  with  the  data  for  each 
modeling  of  a process.  The  number  of  samples  of  X(t)  used  to  obtain 
five  thousand  surges  is  printed  in  the  title  of  the  histogram.  The 
height  of  each  rectangle  of  the  histogram  is  normalized  so  that  the 
summation  of  the  areas  of  all  of  the  rectangles  equals  one. 
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Estimating  the  Autocorrelation  Func-"  ion 
The  estimate  of  the  process  autocorrelation  function  is  also 
plotted  with  the  data  for  each  modeling  of  a process.  ’ To  avoid 
having  to  store  all  of  the  samples  of  the  process,  the  estimate  of 


the  autocorrelation  function  is  based  only  on  the  first  five 


thousand  samples  of  the  process. 


The  estimate  of  the  autocorrelation  function  R(t)  was  obtained 
using  the  biased  estimator 


, N-t 
No  K=1 


where 


N 

'2  1 2 

° ^ ^ K 

N K=1  ^ 


This  estimator  has  the  desirable  property  that  the  variance  of 
the  estimator  does  not  depend  on  t.  Because  R(t)  was  only  calculated 
for  twenty  five  lags  and  the  estimate  was  based  on  five  thousand 
samples,  the  bias  in  the  estimator  is  not  serious. 

The  Surge  Probability  Density  Function 

For  every  process  which  was  simulated,  the  actual  distribution 
of  the  svirge  lengths  p^(ot)  was  approximated  by  a histogram  of  the 
simulation  data. 

Because  the  process  X(t)  was  simulated  by  a discrete  process, 
only  surges  of  integer  length  could  be  detected.  The  rectangles 
are  all  centered  over  integers  and  the  height  of  the  rectangle  over 
integer  I is  equal  to  the  niunber  of  surges  of  length  I,  divided  by 
5000.  Thus  the  width  of  all  rectangles  is  one  and  the  sum  of  all 
rectangle  heights  equals  unity. 

Confidence  Intervals 

The  histogram  of  surge  lengths  represents  an  estimate  of  the 
true  surge  length  distribution.  Because  an  evaluation  of  another 
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estimate  (Denisenko's  p^{ot))  is  made  by  comparing  it  to  the  histo- 
gram, it  is  desirable  to  have  a measinre  of  ho"  accurately  the 
histogram  represents  the  true  surge  distribution. 

To  accomplish  this,  the  5000  sxirges  represented  by  each 
histogram  were  divided  into  twenty  separate  epochs  of  250  surges 
each.  By  examining  the  variation  from  epoch  to  epoch  of  each 
rectangle,  a 95^  confidence  intervsil  was  calculated  for  each 
rectangle.  The  confidence  interval  was  indicated  by  tick  marks  of 
length  one  third  above  ajid  below  the  top  of  each  rectangle. 

In  the  calculation  of  the  confidence  intervals,  it  was  assumed 
that  the  heights  of  the  twenty  rectangles  were  normally  distributed 
about  the  true  value , and  that  the  epochs  were  independent . The 
first  assximption  can  be  partially  Justified  by  saying  that  the 
heights  represent  the  simis  of  repeated  independent  trials.  In  order 
to  have  epochs  with  low  correlation  (if  not  independence),  the 
process  was  run  for  one  hundred  samples  after  the  end  of  each  epoch 
before  starting  to  take  surge  data  for  the  next  epoch.  During  this 
one  hundred  sample  period  between  epochs  data  was  still  taken  for 
the  computation  of  -R"(0)  (to  be  discussed  later),  the  histogram 
of  the  process  probability  density  (not  to  be  confused  with  the 
surge  density),  and  the  autocorrelation  estimate. 

Because  the  process  X(t)  is  modeled  by  the  discrete  process 
Xj^,  certain  quantization  errors  will  be  present  in  the  surge  length 
histogram.  These  errors  are  discussed  in  Appendix  F.  These  errors 
make  it  incorrect  to  say  that,  "the  true  distribution  lies  within 
the  confidence  interval  with  probability  95/?"*  which  would  be  true 
if  the  quantization  problem  were  not  present.  The  confidence 
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interval  is  still  useful,  though,  as  an  indication  of  the  variance 
of  the  histogram.  If  the  confidence  interval  for  a given  rectangle 
is  small  then  the  height  of  the  rectangle  may  well  be  close  to  the 
true  value  of  the  distribution,  and,  at  worst,  it  can  be  said  that 
any  errors  made  were  committed  in  a consistent  manner. 

Epochs  Beginning  in  Mid-Surge 

Breaking  up  the  simulation  into  epochs  introduces  the  potential 
for  another  type  of  error  not  considered  in  Appendix  F.  If  the 
process  is  above  the  siirge  level  when  the  epoch  begins,  then 
the  first  part  of  that  surge  will  have  been  ignored  and  the  resulting 
surge  distribution  will  have  bias. 

To  avoid  this , the  program  checks  at  the  beginning  of  each 
epoch  and  if  it  is  greater  than  x^  then  the  remainder  of  this  surge 
is  ignored  and  the  epoch  effectively  begins  with  the  first  sample  of 
the  next  surge. 

Finding  the  -R"(0) 

To  calculate  Denisenko's  estimate  of  the  s^irge  length  distri- 
bution from  equations  (ll)  and  (12)  or  (17),  the  value  of  -R"(0) 
for  the  process  Xj^  is  required. 

The  autocorrelation  function  of  a second  order  autoregressive 
process  is  known  exactly.^  Let  and  A^  be  the  coefficients  of 
the  second  order  autoregressive  process.  If  the  roots  of  the 

^Gwilyn  M.  Jenkins  and  Donald  G.  Watts,  Spectral  Analysis  and 
its  Applications.  San  Francisco,  Holden-Day,  I968,  p.  166. 
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equation  y - A^y  - = 0 are  and  itg  and  if  these  roots  eire 

real  then 


R(t)  = 


1^1(1  - ' ■ ’^2^^  ” 

2 

(iTl  - 


(20) 


and  if  the  roots  are  complex 

^cos(2'nf  T - ({>  ) 

R{t)  = 2 ^ ^ 

a COS41 

Jv  O 


(21) 


where 


f 


o 


2k 


arc  cos 


and 


1 — 

(t  = arc  tan( 7;  tan  2Kt  ) 

° 1 + 


In  theory  then,  all  that  needs  to  he  done  is  evaluating  the 
second  derivative  with  respect  to  x of  the  proper  equation  ((20)or(2l)) 
aud  evaluating  it  at  zero.  The  problem  is  that  the  function  R(t) 
often  has  a sharp  peak  at  x = 0 and  thus  R(x)  has  no  first 
derivative  at  zero,  and  therefore  the  second  derivative  does  not 
exist  at  this  point.  This  problem  can  sometimes  be  avoided  by 
defining  a new  R(x)  function  in  the  neighborhood  of  zero  with  a 
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"rounded"  peak,  but  the  resxilts  become  very  ambiguous.** 

In  the  continuous  case,  the  quantity  -R"(0)  is  the  expected 
value  of  the  square  of  the  first  derivative  of  the  process.  The 
first  derivative  of  a continuous  process  corresponds  to  the  first 
difference  of  a discrete  process.  This  gives  the  estimator  for 
-R"(0) 


For  N on  the  order  of  fifty  thousand  samples  (which  was 
usually  the  case)  this  estimator  is  accurate.  To  avoid  the  problem, 
outlined  above,  of  derivatives  not  existing,  this  estimator  was 
used  for  the  value  of  -R"{0)  in  the  calculation  of  Denisenko's 
estimate  p^(a). 

Plotting  Denisenko's  Estimate 

To  facilitate  easy  comparison  of  Denisenko's  estimate, 
p^(a),  with  the  histogram  of  the  experimental  data  (a  curve  as  close 
to  the  true  p^(a)  as  is  possible  to  obtain)  the  estimate  p^(qi)  is 
output  as  a smooth  cvirve  on  the  same  plot  as  the  histogram. 

Because  the  process  X(t)  is  modeled  by  the  discrete  process  Xj^, 
it  is  impossible  to  acc\irately  detect  stirges  of  length  less  than 
one.  The  histogram  rectangle  corresponding  to  s\irges  of  length  one 
extends  from  one  half  to  one  and  one  half.  Thus  the  histogram  has 


‘*Petr  Beckmann,  Probability  in  Communication  Engineering. 
New  York,  Harcort,  Brace  & World  Inc.,  19^7,  pp.  22U-226. 
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all  of  its  area  in  the  region  a greater  than  one  half.  The 
estimate  p^(o)  is  of  the  surge  length  distribution  of  the  continuous 
process  X(t)  and  thus  it  has  all  of  its  area  in  the  region  a greater 

A 

than  zero.  If  a direct  comparison  is  made  between  p^(ot)  and  the 
histogram,  one  would  not  expect  a good  fit  of  the  two  curves  even 
if  both  were  error  free  representations  of  the  true  distribution. 

This  is  because  they  have  differing  lower  bounds  on  their  definitions. 

To  avoid  this  problem,  suid  allow  a more  accurate  evaluation  of 
Denisenko's  estimate,  the  smooth  curve  which  is  actually  plotted 
is  a scaled  version  of  Denisenko's  estimate.  Let  p*(o)  be  this 
scaled  estimate.  Then 


p»(a)  = Kp^(o) 


where  K is  chosen  so  that 


fk 


i*(a)  da  = 1 


Thus 


K = 


1 - j p^(a)  da 


and  the  curve  p*(a)  is  only  plotted  for  values  of  a greater  than 
one  half. 


The  Siirge  Cumulative  Density  Function 
In  addition  to  the  plot  of  the  s\irge  probability  density 
function,  a plot  of  the  surge  cxanulative  density  function  is  given 
with  the  data  for  each  simulation  of  a process.  The  experimental 
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data  is  plotted  as  dots  and  Denisenko’s  estimate  of  the  cumulative 
density  is  plotted  as  a smooth  curve.  As  with  the  PDF  estimate, 
the  CDF  estimate  is  scaled  and  only  defined  per  values  greater 
than  one  half. 


CHAPTER  IV 


DISCUSSION 

The  data  from  eleven  different  simulations  is  included  in  this 
chapter.  All  of  the  data  for  each  simulation  is  presented  in  a 
single  figure.'  The  equation  in  the  upper  left  part  of  each  figure 
describes  the  marginal  distribution  of  the  process  Xj^.  A^  and  Ag 
are  the  coefficients  of  the  second  order  autoregressive  process 
(see  equation  (l8))  which  was  used  to  generate  the  process  Xj^. 
Because  the  graphs  are  reproduced  in  small  scale,  the  values  printed 
along  the  axis  are  difficult  to  read.  For  this  reason,  the  extreme 
values  for  each  axis  are  typed  at  both  ends.  Because  all  scales 
are  linear,  the  reader  can  easily  interpolate  between  the  extreme 
values . 

What  constitutes  a "good"  fit  is  a question  which  is  left 
entirely  up  to  the  eye  of  the  reader  and  no  attempt  has  been  made 
to  define  an  index  of  performance  such  as  mean  squeire  error  or  mean 
absolute  error. 

Processes  With  Exponential  Correlation  Functions 

If  the  coefficient  Ag  of  the  second  order  autoregressive 
process  is  equal  to  zero  and  the  process  is  Gaussian,  then  the 
correlation  function  of  the  process  is  exponential.  Figures  3 and  U 
show  that  the  correlation  function  is  also  approximately  exponential 
if  the  process  is  exponential. 


26 


As  discussed  at  the  end  of  Chapter  2,  good  performance  of  the 
estimator  should  be  expected  for  processes  with  exponential  correla- 
tion functions.  The  fit  of  the  estimate  to  the  data  (looking  at 
probability  density  function  plots)  in  figures  1 through  U is  much 
better  for  Gaussian  processes  than  exponential  ones.  Comparison 
of  figures  1 and  2 as  well  as  3 and  U shows  that  the  performance 
remains  roughly  the  same  as  the  surge  level  is  changed. 

Processes  With  Periodic  Components 

The  processes  represented  in  figures  5,  6 and  T cause  the 
estimator  performance  to  be  worse  than  it  was  in  the  exponential 
case.  Again  the  fit  is  better  for  Gaussian  processes  than  exponen- 
tial. Comparison  of  figures  5 and  6 shows  again  that  the  performance 
of  the  estimator  is  not  strongly  affected  by  changes  in  the  surge 
level . 

The  processes  of  figures  8 and  9 have  correlation  functions 
which  are  very  unexponential  in  character.  With  these  processes, 
the  estimator  does  not  fit  the  data  well  at  all  and  does  not  even 
have  the  correct  shape.  These  processes  have  a band  pass  spectrum 
and  so,  while  still  random  and  stationary,  they  have  a definite 
periodic  component.  If  a process  has  a periodic  component  and  the 
surge  level  is  less  than  the  magnitude  of  this  component,  then 
obviously  the  most  probable  surge  length  will  not  be  zero.  This 
intuitive  notion  is  verified  by  figure  8. 

The  form  of  Denisenko's  estimator  guarantees  that  the  estimate 
of  the  surge  duration  distribution  will  be  exponential,  and  thus  it 
will  always  predict  that  the  most  likely  surge  length  is  zero.  With 
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this  in  mind,  it  is  not  surprising  that  the  estimator  does  not  do  a 
good  Job  with  processes  having  periodic  components. 

If  the  process  has  a periodic  component,  and  the  surge  level 
is  adjusted  so  that  the  most  likely  surge  length  is  one,  then  the 
flaws  in  the  estimator  will  occur  for  surge  lengths  less  than  one 
and  the  quantization  discussed  earlier  will  hide  these  flaws  and 
give  a fit  which  is  apparently  very  good.  This  is  what  is  happening 
with  the  processes  in  figures  10  and  11.  These  figures  also  show 
that  the  estimator  performance  is  dependent  on  the  criterion  used  to 
evaluate  performance.  If  for  a certain  application  one  is  concerned 
with  getting  a good  fit  to  the  c\jmulative  density  function  of  the 
distribution  then  the  estimator  would  be  said  to  perform  better  with 
the  exponential  process.  If  instead,  one  were  interested  in  a good 
fit  to  the  probability  density  fiuiction  of  the  distribution  then 
the  estimator  would  be  said  to  perform  best  with  the  Gaussian 
process. 

Conclusions 

It  has  been  seen  that  the  performance  of  Denisenko's  estimator 
of  the  surge  length  distribution  depends  a great  deal  on  the  cor- 
relation function  of  the  process.  In  general,  the  estimator 
improves  as  the  correlation  function  of  the  process  approaches  an 
exponential  function. 

It  is  also  concluded  that  the  estimator  does  not  work  as  well, 
when  applied  to  non-Gaussian  processes.  It  should  be  remembered, 
though,  that  Rice's  estimate  does  not  work  at  all  for  non-Gaussian 


processes,  so  in  this  case  the  extended  Denisenko  estimator  may  well 
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be  as  good  as  one  can  do. 

So  finally,  if  an  application  requires  an  estimate  of  the 
surge  distribution  of  a normal  process,  and  the  correlation  function 
is  similar  to  an  exponential  function,  then  Denisenko's  estimate  may 
well  be  good  enough,  and  is  far  simpler  to  obtain  than  Rice's.  If 
the  correlation  function  is  not  close  to  exponential  and  good 
accxiracy  of  the  estimator  is  required,  then  Rice's  method,  first 
introduced  in  1958 » still  stands  as  the  best  solution. 
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APPENDIX  A 


A COUNTER-EXAMPLE  OF  THE  CLAIM  MADE  IN  EQUATION  (8) 

Consider  the  process  X(t)  which  is  normally  distributed  with 
zero  mean  and  unit  variance  for  all  t.  Let  the  autocorrelation 
series  of  X(t)  be  such  that 

R(0)  = 1 
R(l)  = 0 

and 

R(2)  = r 

euid  consider  surges  of  X(t)  above  level  = 0.  The  claim  of 
equation  (8)  is: 

p(x(t)  > olx(t-i)  > 0)  ^ p(x(t)  > o|x(t-i)  > 0,  X(t-2)  > 0) 

(22) 

Let  X(t)  = X(t-l)  = Xg,  and  X(t-2)  = x^.  Now  because 
COVCx^jXg)  = 0 and  x^  and  x^  are  both  normally  distributed,  x^  is 
independent  of  x^.  Thus 

P{x^  > 0|x2  > 0)  = P(x^  > 0) 

= 1 - <t(0) 


.5 


*♦3 


Now 

P(x  > 0,  X > 0,  X > 0) 

P(x^  > Ojx^  > 0,  x^  > 0)  = p(x2  > 0,  x^  > 0) 

But  again  because  R(l)  = 0 and  Xg  and  x^  are  both  normal,  Xg  and 
x^  are  independent. 


So 


P(x2  > 0,  x^  > O)  = P(x2  > O)  • P(x2  > 0) 


= [1  - '»(0)][1  - *(0)] 

1 1 
■2*2 


= .25 

Equation  (22)  now  says 

P(x^  > 0,  x^  > 0,  x^  > 0) 


.5  > 


.25 


uu 


The  nailtivariate  normal  density  with  zero  mean  is 


P{X)  = 


exp[-  I 2^^  X] 


(25) 


where  d is  the  dimension  and  Z is  the  covariance  matrix. 


In  this  case,  d=3  and 


So 

111  = 1-r^ 


where  A = 


and  B = -rA. 


Now  equation  (25)  gives 


P(X) 


^A  0 bV 

/ 

/M 

exp  - „ (x  ,x  ,x  ) 

(2  ^ ^ ' 

0 1 ■ 0 

'yB  0 A/ 

*2 

1*5 


Evaluation  of  the  product  in  the  argument  of  the  exponential  gives 


i(Ax^  + 2Bx^X2  + Ax^  + x^) 


Returning  to  scalar  notation  and  substituting  into  equation  (2l*) 
yields 


P(x^  > 0,  x^  > 0,  x^  > 0) 


1 

/^_\3/2 


r 2 2 2 

Ax  +2Bx  X +Ax  +x 

^exp 2 '^1  *^2  ^3 


(27t)"'"  /l-r‘ 


2 2 

,«ox>o  Ax  +2Bx^x  +Ax  f<”  X 

— I ^ I i 1 J 

/l-r^  °° 


Using 


4-2 

exp(-  ^)dt  = — 

Jo  Jo 


to  evaluate  the  inner  integral  yields 


P(x^  > 0,  x^  > 0,  x^  > 0) 


r , *^i,r  r 

exp( ~)l  exp  [- 


2BXt  +x_+Ax- 

- '-a-  ^3  <^1 


completing  the  square  in  the  argument  of  the  second  exponential 
yields 


P(x^  > 0,  x^  > 0,  x^  > 0) 


Bx.  2 

Tk 


2 


where 


K,  = 


The  change  of  variable  in  the  inner  integral 


Bx 

z = x_  + 

3 


gives 


P(x^  > 0,  Xg  > 0,  x^  > 0) 


^1  f” 
/a  ■ 


exp 


(A  - 2i)x" 

/I  1 


exp(-  I”)  dz  dx^ 


Bx, 


/A 


Using  the  fact  that 


rexp(- 


^)dt  = — (l  - erf(— )) 


dx. 


gives 


P(x^  > 0,  > 0,  x^  > 0) 


2(2tt) 


r(l  - erf(-^x  ) exp[- 
Jo  ^ ^ 


a/a  - B 2,  , 
X J dx 

2/A  ^ ^ 


and  the  change  of  variable 


/2A 


yields 

P(x^  > 0,  x^  > 0,  x^  > 0) 


P(l  - erf(Bz))  exp(-  (A^  - B^  A)z^)  dz 


which  is  this  standard  form  from  a table  of  integrals^ 


(l  - erf(3x))  exp()j^x^)  dx  = 


where 


r 

J n 


^I.  S.  Gradshteyn  and  I.  M.  Ryzhik.  Table  of  Integrals,  Series, 
and  Products.  Ijth  ed. , New  York,  Academic  Press,  1965,  p.  ^^9  and 
p.  lO^tl.  The  pertinent  equations  are  6.286(l)  and  9.121(7)  • 
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So  finally 

P(*l  ^ -*2  ^ ^3  ^ ” 8BitZ 

where 

Z 

1-z 

Evaluating  this  expression  for  r in  the  autocorrelation 
sequence  equal  to  .36  gives ^ 

P(x^  > 0,  > 0,  x^  > 0)  = .169^+9 

and  recalling  equation  (23)  gives  the  contradiction 

.125  1 .1691J9 

The  value  .169  is  not  surprising  because  .125  is  the  probability 
of  all  three  variates  being  greater  than  zero  if  they  are  all 
independent,  but  in  this  example,  x^  and  are  positively 
correlated. 

In  order  for  this  to  be  a valid  counter-example,  a process  with 
autocorrelation  sequence  R(0)  = 1,  R(l)  = 0,  and  R(2)  = r must  have 
an  extended  autocorrelation  sequence  which  results  in  a positive 
definite  matrix.  As  proved  by  Burg,  the  only  requirement  for  this  to 
be  so  is  that  the  three  by  three  matrix  with  the  above  values  be 


^The  integral  was  first  evaluated  by  a computer  program  using 
trapazoid  rule  integration  from  zero  to  twenty.  This  gave  the 
result  P(x^  > 0,  > 0,  x^  > O)  = .155j  and  the  above  calculations 

are  a verification  of  this. 


r 


I 

» 

! 

1 i 

i t 

I 

! I 
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positive  definite. **  For  r less  than  one,  the  three  hy  three  matrix 
is  positive  definite  so  this  is  a valid  counter-example. 


I 


j 

i 


“^John  Parker  Burg.  Maximum  Entropy  Spectral  Analysis.  A 
Dissertation  Submitted  to  the  Department  of  Geophysics  and  the 
Committee  of  Graduate  Studies  of  Stanford  University  in  Partial 
Fulfillment  of  the  Requirements  for  the  Degree  of  Doctor  of 
Philosophy,  May  19T5,  p-  22. 


L. 
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APPENDIX  B 


NORMALIZATION  OF  THE  AUTOREGRESSIVE  PROCESS 

Let  Yj^  be  the  unnormalized  output  of  the  second  order  auto- 
regressive process. 

where  Z,  are  successive  samples  drawn  independetly  from  a normal 
distribution  with  zero  mean  and  unit  variance.  is  independent 
of  Xj  for  all  k and  j.  The  variance  of  Yj^  is  a function  of  Aj^  and 
Ag  and  will  not,  in  general  be  \inity.  Thus  if  Xj^  = Yj^  the  process 
X^  will  not  be  stationary.  The  task  is  to  find  a normalization 
constant  C such  that  when  Xj^  = CYj^,  Xj^  has  \init  variance  for  all  k. 

\ = * ^2^-2 

2 

and  C will  be  found  such  that  E[Xj^  ] = 1.  Expanding  the  square  of 
Xj^  gives 


* “A-2*k  * 


Collecting  terms  and  evaluating  the  expected  values  gives 


where 


' * 'VIVA-2’ 


Using  the  stationarity  of  Xj^  gives 


CK 


1-CA^ 


Substituting  this  result  into  equation  (27)  yields 


o o o 2A-  A_C 


or 


1 - CAg  = C^aJ(1-CA2)  + Agd-CAg)  + 2A^A2C  + 1 - CA2 


Regrouping  like  powers  of  C gives 


0 = C^(A^A2  - A2  - A2)  + C^(A^  + A2  + 1)  + CAg  - 1 
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(28) 


The  process  generating  subroutine  solves  this  cubic  equation 
for  the  three  roots  of  C,  picks  the  smallest  positive  real  root 
and  uses  this  value  to  normalize  the  process. 


APPENDIX  C 


INITIAL  VALUES  OF  THE  AUTOREGRESSIVE  PROCESS 

Often  when  generating  an  autoregressive  process,  one  sets  the 
initial  values  equal  to  zero.  This  gives  a process  which  is  not 
in  steady  state  from  the  beginning,  but  has  a start  up  transient. 
Instead,  the  initial  values  X*  ^ and  X^  ^ 'the  second  order 
autoregressive  process  should  be  chosen  according  to  the  formulas 


^-1  = 


(29) 


and 


x»_2  = RZ^Wi-r  Z, 


(30) 


where  the  values  and  Z^  are  independent  samples  from  a normal 
distribution  with  zero  mean  and  unit  variance.  R is  the  covariance 
of  Xj^  ^ and  Xj^  ^ which  was  shown  in  equation  (28)  to  be 


KA, 


R = 


1-KA, 


By  inspection,  X*  ^ and  X|^  ^ have  zero  mean,  and  X*  ^ has  unit 
variance.  The  variance  of  X*  ^ -s 
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Var(X*.2)  = 

= E[(RZ^  + Z^)^]  - (E[RZ^  + Z^)^ 

= E[R^Z^  + aRZ^/^Zg  + (1-R^)Z2]  - (RE[Z^]  + /C?e[Z2])^ 
= R^E[Z^]  + 2R  A-^  E[Z^Z2]  + (1-R^)E[Z2] 


= R^  + (1-R^) 
= 1 


So  the  variance  of  2 is  correct. 
The  covariance  of  X^  ^ and  X^  ^ 


C0V(X-_i,X*_2)  = E[X*_^,X*_2]  - E[X»_,]e[X*_2] 


= E[Z^(RZ^ 


+ /l-R^'Zg)] 


= E[RZ^  + /i-R^'z2Z^] 

= RE[zJ]  + YXZ^^ 

= R 


Thus  it  has  been  shown  that  if  X*  ^ and  X*  ^ are  picked 
according  to  equations  (29)  and  (30)  then  they  have  the  proper 


means,  the  proper  variances,  and  the  correct  covariance. 


APPENDIX  D 


COMPUTATION  OF  THE  MAPPING  FUNCTION  G(x) 


Equation  (l9)  gives  the  napping  of  a normal  process  into  an 
exponential  process.  The  computation  of  this  function  was  required 
every  time  that  an  exponential  variate  was  needed.  The  program 
generated  about  50,000  exponential  variates  for  every  simulation, 
so  an  efficient  way  of  computing  the  error  function  was  required. 
The  error  function  was  approximated  in  the  following  way.^ 

2 

erf(x)  = 1 - (a^t  + a^t^  + a^t^  + aj^t^  + a^t^)e”* 

where 


p = .3275911 
a^  = .25U829592 

ag  =-.28UU96736 

a^  = l.U2l4l37l+l 

aj^  =-1.U53152027 

a^  = 1.06lit05i*29 

and  the  magnitude  of  the  error  e(x)  is  less  than  1.5’10  * 


^Milton  Abramowitz  and  Irene  Stegun.  Handbook  of  Mathematical 
Functions , National  Bureau  of  Standards,  I96U,  p.  299* 


I 
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APPENDIX  E 


THE  MEAN  OF  DENISENKO'S  ESTIMATE 


Equation  (ik)  gives  Denisenko's  approximation  to  the  surge 
length  distribution  as 


p (a)  = Ae 

T 


-Aa 


where 


A(x^)  = 


/-R"(0) 


(x  ~ \i] 


v2 


u - vT  '=‘1’ ' rs- 


2 [1  - 4> — ] 


2a 


and 


4>(z) 


z 


dt 


So 


E[t]  = T- 


2Tr[l  - $ { 


° ■)] 


/-R"(0) 


a - K 

exp  [+  2 

2a 


exp  ] {2[1  - $ (-2_)]}  (31) 


/-R"(o")  2a^ 


With  some  algebra,  the  function  $ can  be  expressed  in  'terms  of  the 
error  function: 


I 


substituting  this  into  (3l)  yields 


(x  - y)  11  *0  ” ^ 

Efn  (a)l  = exp  ] {2[l  - ^ - erf  ( — )] 

^ /IFW  2a^  ^ ^ 


(x.  - y)' 


x_  - y 


Efn  (o)]  = exp  [ — - — ?; ] {1  - erf  : — )}  (32 

^ /-R"(0)  2a  a/2 


And  now  for  comparison,  the  exact  mean,  equation  (2),  is  rewritten 


as  follows 


(x  - y)  X - y 

Etp  (a)]  = - exp  [ — - — 5 ] {erfc  ( — — )} 

/-R"(0)  2a  a/2 


E[p  (a)]  = ■; — exp  [ 5 ] {l  - erf  ( — — )} 

✓-R'HO)  2a"  a/2 


Comparison  of  this  expression  with  equation  (32)  shows  that 
Denisenko's  approximation  gives  an  unbiased  estimate  of  the  mean. 


APPENDIX  F 

ERRORS  DUE  TO  QUANTIZATION 

The  process  X(t)  is  modeled  hy  the  discrete  process 
can  be  thought  of  as  the  sampling  of  X{t)  at  integer  time  increments. 
A possible  segment  of  the  process  X(t)  is  shown  in  figure  12  where 

tick  marks  indicate  sampling  times  and  the  socis  is  the  surge 

threshold.  The  following  list  details  some  of  the  possible  errors 
due  to  quantization.  Numbers  refer  to  the  numbers  on  figure  12. 

1.  A short  surge  is  ignored. 

2.  A siirge  of  length  .95  is  ignored. 

3.  A short  surge  is  recorded  as  having  length  one. 

U.  A surge  of  length  1.9  is  recorded  as  having  length  one. 

5.  A s\irge  of  length  .5  is  ignored. 

6.  Three  separate  surges  of  lengths  2.5,  3.5,  and  2 are 
recorded  as  a single  surge  of  length  nine. 

7.  A series  of  short  surges  is  ignored. 

8.  A series  of  short  surges  is  recorded  as  a single  s\irge 


of  length  four. 


I 


i 


APPENDIX  G 

THE  COMPUTER  PROGRAM 


The  process  modeling,  the  calculations  of  the  estimators, 
and  the  output  of  data  were  all  done  with  a FORTRAN  program  written 
for  the  CDC  "RUN"  compilor.  All  graphical  output  is  done  with  the 
CDC  280  microfilm  plotter,  using  subroutines  in  the  CUGLIB  library. 

In  addition  to  the  four  plots  shown  in  figures  one  through 
eleven,  the  program  outputs  an  additional  plot,  an  example  of  which 
is  given  in  figure  13-  This  plot  is  for  the  data  of  simulation 
number  six.  The  plot  is  the  cumulative  distribution  function 
subject  to  the  mapping 


X*  * — In 
^ A 1-X 


where  A is  given  in  equation  (lO)  of  the  text,  K is  the  normalization 
constant  derived  in  Appendix  B,  X is  a point  of  the  cumulative 
distribution  fiinction  and  X*  is  the  mapped  point  which  is  plotted. 

The  nature  of  this  mapping  is  such  that  Denisenko's  estimate 
becomes  the  straight  line  y = x.  Thus  to  see  how  well  the  experi- 
mental data  fits  the  estimate,  one  looks  at  how  straight  a curve 
throxigh  the  points  is.  The  line  y = x is  also  plotted  for 
reference.  It  should  be  pointed  out  that  no  additional  information 
is  contained  in  these  plots  which  is  no.t  already  in  the  regular 


cumulative  density  function.  For  this  reason,  it  is  not  included 


6] 

with  the  regular  data  in  figures  one  through  eleven,  but  it  does 
give  a quick  indication  of  how  well  the  estimator  performs  for  a 
given  process  and  the  routine  has  been  left  in  the  program. 

Also  included  in  the  program  are  subroutines  which  sex'arate 
on  the  microfilm  the  plots  of  one  simulation  from  those  of  the 
preceding  simulation  as  well  as  that  following. 

Each  simulation  requires  one  data  card  which  contains  the 
surge  level  in  columns  two  through  eleven,  the  coefficients  of 
the  autoregressive  process  (A^  in  columns  twelve  through  twenty-one 
and  Ag  in  col\amns  twenty-two  through  thirty-one)  and  a one  in  column 
one  for  a normal  process  or  a zero  in  column  one  for  an  exponential 
process.  Preceding  these  cards  is  a seed  card  for  the  random 
number  generator  which  is  any  number  in  octal  twenty  format.  The 
last  card  in  the  deck  should  have  a single  nine  in  col'jmn  one. 
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.BESLMIUBlEIcaf 


PROGRAM,  SEP  I ES  ( I NPliT  , OUTPUT  , FlLKPL  i T.'»PE5  = IN?UT  i TAPES-aUTPUTI 
01  MENS  I CN  OMARI  iOO)  , S'JRG  ( 20  i 250  I , OUTPUT  I 5000 ) 

REAOIStll  SEED 
1 FORMAT  (C20I 
a*RANFI SEED) 

DO  2 t*i.iao 

REA0I5.T)  IFvnRM.SLEV,  A1.,A2 
3 FORMAT! U ,3F10.a| 

IFIlFNQPR.cO.'il  OC  Tt*  A 
CALL  CMINiTIO.) 

CALL  INCISP 
CALL  frame 
CALL  ARROmO 

CALL  OHAINl Al,A2,SLuV,CMAR,SURGtOUTPUT, IFNORH) 

2 CONTINUE 
A CONTINUE 
STOP 
ENO 


SUBROUTINE  ARROWUt A1,A2,SLEV, IFNORM) 
CALL  INCISP 

CALL  MAF(0.,l.,0.,l.fO.,l.,0.,l.) 
CALL  SETLINE  U ) 

XSHFT=0. 

5 CONTINUE 
00  1 I3lt3 
Rt>l 

X«. 135-.006*{RI-1. )*XSHFT 
YST0P*.897-.006*IRI-l . ) 

1 CALL  LYNE (X, .05 tX.YSTOP  ) 

00  2 1=1,3 

RI  = I 

X*.lAl*.005*(RI-l. I+XSHFT 
YST0P=.897-.006*(«I-1.  ) 

2 CALL  LYNE 1X,.05,X,YSTCP» 

00  3 1=1,6 

RI*1 

YT0P».9+.006*IRI-l. » 

YBTM=YTCP-.l 38 
XLFT=XSHFT 
XmD=.138»XSHFT 
XRHT=.276*XSHFT 

CALL  LYNEIXRHT, YPTM.Xmic.ytOPI 
3 CALL  LYNE IXM 10, YICP.XLFI , YBTMI 
IF  USHFT.NE.O)  UC  TC  a 
XSHFT=.72A 
GO  TO  5 
A CONTINUE 

CALL  LAEPLT  I Al , A2 , SLEV , IFNCRM I 

RETURN 

End 


y&JIMMJM  » 


SUBROUTINE  ARRCWO 
CALL  INCISP 

CALL  MAPtO.i l.f  O.t l.tO»t 
CALL  SETLINEIU 
XSHFT=0. 

5 CONTINUE 
00  1 I=>lt3 
RI^I 

Xa.l35-.006*(RI-l. » tXSHFT 
YSTGP=.l0  3+.005*(RI-l.  ) 

1 CALL  LYNE ( X, .95 t Xt YSTQF ) 

00  2 I«lt3 

R|«I 

Xa.lAl*.00&* (RI-1.  )+XSHFT 
YSTCP=.l0T».0C6*(Rl-l.> 

2 CALL  LYNcIX, .95tX,YST0P) 

DO  3 1=1 i6 

Rl=l 

YBTP=.l-.006*(R I-l. I 

YrOP=YBIP*.l38 

XLFT=XSHFT 

XMIC=.l3e»XSHFT 

XRHT=.276+XShFT 

CALL  LYNE (XRHTt VTCPf XMIOfYSTHI 

3 CALL  LYNEtXPID.YKTM.XLFrtYTGPI 
IF USHFT.NE.O)  00  TO  4 
XSHFT=.724 

GO  TO  S 

4 CONTINUE 

CALL  OPTION! 1,0, 3iOtOI 

CALL  CSTHING! .3945, .6381 , IZHPLEASE  POTS-J 
CALL  CSTRINGI .3594, .5791, ISHTHE  FCLLOV, ING$.  I 
CALL  CSTRINGI .3594, .5264, 15HFCUR  PLOTS  ONF.) 
CALL  CSTRINGI. 413, .4736, IChTHE  SAPEt.) 

CALL  CSTR  1NG( .3828, .4209, 13H3  X 10  INCH*. I 

CALL  CSTRINGI .45 3 I , . 3682 , 7HPR INTt. I 

CALL  frame 

RETURN 

END 


SUBROUTINE  LABPLTI A1,42,slEV, IFNCRMI 
CALL  NAP(0.,l.,0.,2.,0.,l.,0.f 1.) 

CALL  OPTICN  11,0,3,0,0) 

CALL  CSTRING(.383,.6,6HAl  =$.) 

CALL  CSTRI\GI.383,.5,6HA2  =1.) 

CALL  CSTRINGt-453,.35,3h=$.) 

CALL  CNUM6R  ( .5,  .6  , A 1 ,4|-F8.4  ) 

CALL  CNLPRRI .5,  .5,A2,4hF8.4) 

CALL  CNUMBR(.5,  .35,SLEV,4FF8.4) 

CALL  CSTR  ING ( .3  24, . 37 3, 6FSURG S. ) 

CALL  CSTRING(.312,.326, THLEVELi. ) 

IF(  IFNCR''.£C.l)  GC  TQ  I 

CALL  CSIRING(.3C3,.  74  7,  I 3h EXPONENT  I ALT. 1 

GO  TO  2 

1 CALL  CSIR  ING ( .44 1 , . 74  7,8HN0RMAL». I 

2 CALL  CSTR ING( .4 3, . 7 ,9hPRCCESS$. ) 

CALL  FRAME 

RETURN 

ENO 


i 
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SUBROUTINE  CVAIN  I A I , A2  , SL  EV , CM  AR  , SLRG,  OUTPUT  , IFriOR^fl 

UIPENSICN  Cva«(  lOO»,SURG<iOf250),OUrPUT(50CCI 

iNOCUTaO 

ITCTL=5000 

IRUN=TTCTL/20 

RUNalRUN 

Ra.95 

BETAaS. 

ALPHAaO. 

ISRGaO 

ISFSTaO 

llRUNal 

IPOCH=l 

ICNTaO 

IRUNCaO 

RTCNTaO. 

RTAUaO. 

0LST=0. 

Clal. 

DO  2 Ial,100 
OHARdlaO. 

DO  2 Jal,20 

2 SURS(J.I)aO. 

3 continue 

CALL  process  (Cl, alpha, beta, out,  IFN0RH,Al,A2l 
RTCNT=RTCNT*1. 

RTAUaRTAU+tCLST-CUT )*(CLST-OUTI 

OLSTaOUT 

INCCUT=INDCUT*1 

IFtINOQUT.GT.5000)  GO  TO  11 

OUTPUT t INOCUT)=CUT 

60  TO  12 

11  INDCUT=5000 

12  IFt  IFNORv.ec.O)  go  TO  13 
I0UT=OUT*5.*50.5 
IF(IOUT.LE.l)  ICUT=1 

GO  TO  14 

13  I0UT=(CUT*2. )+l. 

14  IF( lOUT.GE.lOO)  ICUT=100 
DMAKdOUT  laDMAR  ( ICUD  + l. 

IF tCUT.GE.SLEV)  ISRG=1 
IFdSRG.EC.l.ANC.  IlRUr+.EC.l  I ISFSTal 
llRUNaO 

IFICUT.LT.SLEV.ANC.ISRG.EC.li  GO  TO  4 
IFdSRG.EO.OT  GC  TO  3 
ICNT=ICKT*1 
GO  TO  3 

4 IFdSFST.6C.II  GC  TO  17 
IFdCNT.GE.lOO)  ICNTalOO 
SURGdPCCH,  ICNT  laSURGt  IPCCH,  ICNT)*1. 

IRUNCal aUNC*l 

17  ISFSTaO 

ISRGaQ 
ICNTaO 
GO  TO  3 

5  IPCCHaiPCCH^l 

IFdP0CH.GE.21l  GO  TO  IB 
DO  19  Ial,ico 

CALL  PROCESS  t C I , ALPHA, BET A, OUT , IFNCRM, A 1 , A2 1 
RTCNTaRTCNT* 1. 

RIAU*RTAU+tCLST-CUT )*(OLST-OUT) 

OLSTaOUT 
iNliOUTa  INDCUT*! 

IF  dNDCUT.GT.5C00)  GC  TC  20 


idi 
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OUTPUT! I\UUUTI=0UT 
GO  ro  29 

20  INOCUT*5000 

29  IF ( IFNQRP.EC.O)  GO  TO  27 
I0Ur»0UT*5.+50.5 
IFtlOUT.LE.U  ICUT=1 
GO  TO  28 

27  10UT=(QUT*2.  I*!. 

28  JFlIOUT.CE.lOO)  ICUT=100 
OMARIICUT)=OMAR(ICUT)*1. 

19  CONTINUE 
llRU.N>t 
IRUNC^O 
GO  TO  3 

18  00  21  (>ltlOO 

SUI**0. 

00  22  J=li20 

22  SUP=SUK*SURG{Jt n 
XBA9>SUP/20. 

IFIXBAR.EC.O.I  GO  TO  25 
ci*0. 

DO  23  J*1.20 

23  CI«ISURGIJt  I )-XlsAR)*I  SURGI  J,  I )-XBAR» 

Cl-CI/19. 

CI>.46AA*SCRT(Cn 
GO  TO  26 

25  SURG(2,M’°0. 

GO  TO  21 

26  SUR3(2t  n=CI/RUN 

21  SuRG(l,n=X8AR/RUN 
HRITE<6.8)  RTCNT 

8 F0RPAT(*1  THE  MARGINAL  CISTRIBUTICN  OF  THE  PROCESS  BASED  ON  •rFB 
1.0. ♦ SAPPLES*) 

KRlrE(6.ISI  A1,A2 

15  FOR>>AT(*OA1  = «,F8.A,*  A2=*.F8.4) 

WRirE(6.16) 

16  FORPAr(*0*) 

DO  6 IslflOO 

IF ( IFNORP.cC.l)  CHARI n = 5. ♦CM AR{  I) /RTCNT 
IFI  IFNORH.EC.O)  CHARI  I I = 2. ♦CHAR ( I ) /RTCNT 
6 WRITE(6.7>  I.CHARII) 

7 FQRHATI I5,F20.8 I 
RTAU*KTAU/RTCNT 
SURGI3.IJ=irCTL 
SURG(3,2  »=RTCriT 
SURC(3,3) =SLEV 
SURGO.AI^AI 
SURGI3r5)>A2 
SURGt3.6)=INCCUT 
SURGI3.7)=lFNORM  . 

CALL  OUTPLCT/DHAR,SURG,OUTPUT,IFNORH,RT4U) 

RETURN 

END 
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SUBROUTINE  CUTFLCr  ( CHAR,  SURG  tOUTFlj  T , IFNOTM  .RUUI 
CIKIENSICN  CMARI  l00)tSU«G«20t250)fUUTPUT(5000) 

CALL  SETLINEdI 
CALL  INOISP 
61 0^0* 

00  SA  1^1,100 

54  IFICMARIII.GE.eiG)  BlC^OyARt  t) 

VMAX*S. *810/4. 

11*100 

55  IFlOMARIin.GT.O.I  GO  TO  56 
11*11-1 

GO  TO  55 

56  R11*II 
IBIG*II 
XKAX*RI 1/2. 

IF(IFNCRK.EQ.l)  XK AX=R 1 1/5 .-9 .9 
XNIN*0. 

ISFL=l 

IFI  IFNORM.EC.O)  GC  TO  60 
11*1 

62  IFICKARl Ilt.GT.O.)  GO  TO  61 
11*11*1 

GO  TO  62 
61  R1I*1I 
ISKL*II 

XMIN=tRII-l. 1/5. -9. 9 
60  CALL  OPTION! 1.0,2, 0,01 
1*SURG(3,2> 

CALL  CSTRING(242,973,4lhPARGIMAL  OISTRIBUTICH  OF 
CALL  CNUM8R  (650. 97  3,  U2FI 6 ) 

CALL  CSTRING(500,75,3HX*.> 

call  MAPS(X?'IN,Xa<AX,0.,YPAX,  .15,.95,  .15,  .91 

IFI  IFNOa.r.EC.O)  GO  TO  63 

YTCP=OMAR(l) 

CALL  LYNEIXPIN, YTCP,XMIN,0.) 

63  00  57  I*ISPL,IB1G 
R1*I 

XSTART=(RI-l.)/2. 

XENC*RI/2. 

IF( IFNQRP.EC.l)  XSTART*tRI-l.  )/5.-9.9 
IF! IFMORP.EC.l)  XENC^XSTART *.2 
YTCP=OMAR(I) 

IF(I.EO.IOO)  GC  TO  59 
YNT=OMAS ( I + l ) 

IFIYTOP.GE.Y.NTI  GO  TO  59 
CALL  LYNE(XEND,YTCP,XENO,YNT) 

59  CALL  LYiVeiXSTARr.YTCP.XENC.YTOP) 

57  CALL  LYNe(XEN0,YTCP,X£N0,0.» 

CALL  OPTION! 1.0,2,1,01 

CALL  CSTRING!75,400,l4h0E.NSlTY  OF  XS.) 

CALL  FRANE 
IN0CUT=SURG(3,6» 

A1*SURCI3,4) 

A2*SURG!3,5I 

CALL  SPIM!CUTPUT, INROUT, A1,A2) 

CALL  INCISP 

eiG*o. 

00  50  1*1,100 

so  IFISURGIl.n.GE.aiG)  BtG*SURG!  1,  I ) 

YMAX=5. *810/4. 

IF ( IFNORM.EG.OI  GO  TO  64 
X*SURG( 3, 31/ 1.4 14213562 
CALL  ERFNCTN!X,V) 

PHI*.5*Y/2. 


SAMPLES*. I 


j 


) 

i 


i 


J 
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fCTR^JRTAO/  t6.2ima5i*(  1-PHU  H »EXP l-.5*SURC ( 3i  3 ) • SU-tli t 3 1 J1  I 
GO  TO  66 

66  FCTR-SOaTJRr.VJ/3.) 


66  CORl*EXP(FCTR/2.) 

YrOP=FCTR»ex?(-l.*FCTR*.5)PCORT 
IF« YTOP.CE.eiG*  YPAX»5.*YTOP/6. 

11-100 

51  IFISURGll.ni.GT.O.I  GC  TO  52 
II-ll-l 

GO  TO  51 

52  Rll-tl 
XHAX-RI t«.5 
I«SliRG(3.  U 

CALL  OPT10N(l,Cf2.0,0) 

CALL  CSTR  INCt233,973,6lHNnRMAHlEU  DISTRIBUTION  OF  SURGESS.I 

CALL  CNt«;tR(  707,973.  If2HI5) 

CALL  CST«ING(956, 75, lAhSURGE  LcGNTHl.l 
CALL  MAPSI.5,XPAX,0..YMAX, .15,.9S,.15..9) 

YTOP=SUPC(l,l) 

CALL  LYNE(.5,0. ,.5,YT0P) 

00  5)  I =>1.11 
RI-I 

XSTART-BI-.5 
XEND-RI *.5 
YTCP=SURG( 1. n 
IFII.EQ.100>  GO  TO  58 
YNT-SURG(1. I*l) 

IFIYTOP.GE.YNDGC  TO  58 
CALL  LYNEIXENO.YTOP.XENC.YNT) 

58  CALL  LYNE (XSTART.YTCP.XEND.YTOP) 

CALL  LYNEIXENO. YTCP.XENO.O.) 

XSTART=XSrART+l./3. 

XEND-XSTART+1./3. 

YTCP-SURGd.  n + SURG(2,II 

CALL  LYNE (XSTART.YTQP.XENU.YTOP) 

YTCP-SURGd,  1 )-SURG(2, 1 ) 

53  CALL  LYNE(XSTARr,YTQP,XENC.YrOP) 

XlNC-RII/500. 

XSTART-.5 

YTCP«C0Br*FCTR*EXP(-l.*FCTR*.5) 

00  65  1-1.500 

Rl-l 


65 


15 


16 


10 

11 


xeNC-.5*RI*XIN'C 

YNT-FCTP*CCRT*EXP(-1.*FCTR*XEND> 

CALL  LYAEIXSTARr.YTOP.XENC.YNTI 

XSTART-XENO 

YTOP-YNT 

CALL  OPTICN  (1.0. 2.1.0) 

CALL  CSTRING(75.62C. 18KNUPBER  OF  SURGES!. > 

CALL  FRAME 
WRITE(6.9) 

FORMAT!*!  THE  SURGE  DISTRIBUTION  OF  THE  PROCESS*) 
WRirEI6.15)  SURG(3.A),SURG(3.5) 

FURMAT(*0Al-*,F8.4.*  A2>*«F8.6I 

MKITE<6.16) 

F0RMATI*0*) 

00  10  I-l.II 
RI-I 

T1«FCTR*EXPI-1.*FCTR*RI ) 

I2«T1*CCRT 

WRl TE(6.1 I)  I ,SUPG( 1. I ).SURGI 2. I ).T1.T2 

FORMAT! 15  ,F15.8.*  *08-* . F 1 5. 8, • T1-*,F15.8,*  r2-*,F15.81 

SURG!6,1)*FCTR 

SURG!6.2>-CCRT 


i 


i 
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SUBROUTINE  Cu^  ( SURG . 1 1 . ftl AU  I 
OIMCNSICN  SU><G(20.250) 

XNAX>II *l 
CALL  [NCtSP 

CALL  OPI ICN( 1.0,2,0.01 

CALL  CSTRING(<.56,75,  lAHSUKGE  LEGNTH.I 

I-SURG(3,ll 

CALL  CNUwesi 707,973, I, 2H15I 

CALL  CSTRIfiG(208, 77  l.'ilHCUHULAriVE  DISTRIBUTION  OF 
CALL  OPTIOM  t.0, 2.1.0) 

CALL  CSTRING(75,'i2O.ieHNUP0ER  OF  SURGESS.I 
CALL  MAFS(0.,XPAX,0..1.2,.15, .95,. IS,. 9) 

Y»0. 

00  1 I>1,II 
X»1 

Y>iY»SURGn,n 

1 CALL  XSTRINC(X,Y,3H*i.) 
fCTR=SURGK,  1) 

C0RT=SURG{A.2) 

HR1TEI6.6)  RTAU 

6 FORPATI*!  THE  VALUE  OF  -R  DOUBLE  PRIME  10)  °*,F15.8I 
WRITEI6,9)  FCTR 

9 F0RPAT(*0  THE  VALUE  OF  ALPHA  -*,F1S.8I 
NRITE(6,10)  CORT 

10  FORMATING  THE  VALUE  OF  THE  INTEGRAL  TO  ONE  HALF  =*,F15.8I 
RIUII 

XINC»(RII-.5)/500. 

XSTART*.5 

YT0P»0. 

00  2 1> 1,500 
RI«I 

XEND=RIAXINC*.5 

YNT»l.-CCRT*FXP (-l.*FCTR*X£NO  ) 

CALL  LYNEIXSTART,YTCP,XENC,YNT) 

XSTART*XEND 

2 YTCP*YNT 
CALL  FRAME 
AlsSURGtS.A) 

A2xSURGI3,S) 

IFNCRM=SURGI3,7) 

SLEV»SU?C(3,3) 

CALL  ARRGWUI A1.A2,SLEV, IFNORMI 

CALL  INCISP 

CALL  OPTION! 1,0,2.0.01 

CALL  CSTRING(A56,75,19HSLaGE  LECNTHS.) 

I>SURGI3.1) 

CALL  CNUM0R(707,973, t,2HI5 » 

CALL  CSTRINGI233,97),<,lhCUMULATIVE  DISTRIBUTION  OF 
CALL  OPTION!  1,0,2, I ,0) 

CALL  CSTRING(75 .AEO.ISHNUMBER  OF  SURGES!.) 

CALL  MAP(0.,1.,C.,  1.2,  . 13. .95,.1S. .9) 

CALL  LYKE(0..0.,0., 1.2) 

CALL  LYNE(0.,0.,  l.,0.  ) 

CALL  LVNE(0.,0. .1.,  1.) 

CALL  OPTICM 1,0, 0,0,0) 

DO  3 t'1,13 

RI-l-l 

R|>R1/10. 

CALL  LVNE(-.003,RI, .003,RI  ) 

3 CALL  XNUMBR(-.030,R1,RI,AI-F3.1I 
XLSr>0. 

LO  A I«l,  II 
Rl-I 

X«l.-CCRT*£XPI-l.AFCTR*RI) 


SUCCESS.) 


SURGES!.) 
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OIF*X-XLST 

ifix.ce:..?b>  gc  tc  t 
IF (OIF.lt. .oil  GO  TC  A 
XLST«X 

CALL  LYNE(X,-. 003.x. .003) 
XST»X 

IFII.GT.9)  XST=X-.005 
call  XNL'MPR(XST.-.0Z2,(,2H2I 
A CONTINUE 

CALL  OPTICM  1.0.2, 0.0) 

7 V«0. 

OO  9 1=1.11 
ltl«l 

Y»YASURG( 1. I ) 

X«l.-CCRT*ex?(-l.*FCT.^*Rn 

5 CALL  XSTPING  (X.Y.3H»t.) 

CALL  FR4f<E 

RETURN 

ENO 


SUBROUTINE  PROCESS  I C I . ALPHA. BET A. CU T. I FNORH. Al .A2 ) 

IF  ICl.EC.O.)  GC  TO  1 

Vl»A2*tAl»Al-A2»A2-l.) 

V2=Al*Al*A2*A2+l. 

IFIVI.EC.O.)  GO  TC  lA 

P=V2/Vl 

C=A2/V1 

R»-l./Vl 

VA*|3. 1/3. 

VB=(2.*P*P»P-9.*P*Q*27. ARI/a?. 
RAD*VR*Ve/A.+VA»V4*VA/27. 

IF{R40. LT.O. ) GO  70  5 
IFIRAO.EC.O.  I GO  tc  6 
ARG=-Va/2.»S2RT(RA0) 

IF(ARG.LT.0. I CO  TO  7 
ACAP=ARGA*( 1 ./3.  ) 

GO  TO  8 

7 ACAP=-1.aI-I.*ARG)**(1./3.) 

8 ARG*-VB/2.-SCRT (RAD) 

IF  (ARC. LT.O.  ) GO  TO  9 
BCAP:>ARG**II./3.I 

GO  TO  10 

9 eCAP=-l.*{-l.*ARG)**(l./3.) 

10  FACT=ACAP*BCAP 
GO  TO  11 

6 ARG=-VB/2. 

IF (ARG.LT.O.  ) GO  TO  12 
ACAP=ARGAA(l./3.) 

GO  TO  13 

12  ACAP=.-l.*(-l.*ARC5**(l./3.) 

13  FACT*-ACAP 
GO  TO  11 

5 XIP*SORT{-l.*RACI 
PKB=-V0/2. 

RMAG«SCRT(PK8*PKB*XIN*Xl)') 

ANGL-ATAN2(Xlr',PK0) 

REL=>RMAG*»(  1 ./3.  )*CCS(ANGL/3.  ) 

X1FAJ=RPAG*»{ 1./3. )*SIN(  ANCL/3. ) 

FACT1=2.*REL 

FACT2=-FACTl/2.-Xlr«Aj*1.732050a08-P/3. 
FACT3=-FACTI/2. 1 MAj*l . 7 32050808-P/3. 
FACrl»FACTl-P/3. 

FACT»1.E50 

IF (FACT  I .GT.n.. ANC. F act I.LT. FACT)  FACT  = FACT1 
IF IFACT2.GT.0.. ANC.FACra.LT.FACT  ) FACT  = FACT2 
lF(fACr3.Gr.C..ANC.FACr3.LT.FACr)  FACI»FACT3 


BtSLAYA'WLCOPY 

CO  TO  15 

11  FACT»F*CT-P/3, 

CO  TO  IS 

14  FACrxl./$C9r( 1. ♦4l*4l*A2«A21 

15  R1=«FACT»AI/(  l.-FACT*42) 

IF(FACT.GT.l.  ) Vi^lTE(6,2)  Al,A2fFACT 

2 F0RHAT(*0  FOR  A1  = •,F'3-Ct*  ANC  A2  =*tF9.6t*  THE  SMAHJEST  FCSITl 
IVE  ROOT  WAS  *,tl5.3) 

C22*SQRT(  l.-Rl*i<l) 

C1»0. 

NCH«1 

CALL  NCRWIXNCRH.NCHI 
XM2»XNCRH 

CALL  NCRW(XNCRH,NCH» 

XHl*RlA*M2'-C22»Xf«CRM 
1 CALL  NORM(XNOR« ,NCH) 

X>{A1*XN1  *A2^XH2*Xt42liH)*FACT 

XHZ^XHl 

XM1>X 

IF(IFNORW.EQ.O)  GO  TO  3 

CUT>X 

GO  TO  4 

3 X*X/l. 414213562 
CALL  FRFNCTNtX.Y) 

0UT»0ET4*(-1.*ALQGC.5-.5*Y l+ALPHAI 

4 CONTINUE 
RETURN 
ENO 


SUBROUTINE  NCRF  tXNCRM,NCHI 

PI»3. 141592654 

IF  INCH.EU.O)  GC  to  1 

Ul^RANFIOl 

y2*RANF(0» 

XRC-sJSCRT  I -2  .*ALCG(un  ) I*COS(  2.*PI*U2» 
XRE=(SCRT (-2.*AlCG(UI ) ) )*SIN( 2.»PI*U2» 
XNORM«XRO 
NCH=0 
GO  TO  2 
1 XNCRH^XRE 
NCH>1 


2 CONTINUE 
RETURN 
ENO 


SUBROUTINE  ERFNCTN  (X,YI 
P*. 32759 1 1 

f Al«. 254829592 

i A2«-. 294496736 

A3»l. 421413741 
A4*-l.453l52C27 
A5*l. 061405429 
SINE>1. 

IF  lX.GE.O.)  GC  TO  1 
SfNE>-l. 

X*-1.*X 

1 T*l./(l.4P*X) 

POLV*Al  *T4’A2«T*r«Al*T*r*T>A4«r**4.'»A5*T**5. 
V»SINE*( l.-PCLY*tXPt-X*X )) 

X«X«SINE 

RETURN 

ENO 
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SUBROUTINE  SPTf*  I CU IPUT , INOOUT , A I,  A2  J 

OIKENSICN  CUTouKSOOOl 

SUP«0. 

W(IITE(6.A)  INDC'J’’ 

4 FORPATC*!  IHt  CORRELAriON  FUNCTION  OF  THc  PROCESS  BASED  ON  *,15 
1,«  SAMPLES*! 

MRl1E(6,dl  Al,A2 

8 FORMAFI «OAl=*,Fa.4, ♦ A2>**F8.4) 

URITE(6t6l 

6 FORMAT(*OLAGS  CORRELATION*! 

RNOOUT»INCCUr-l 

CALL  INCISP 

CALL  OPTiCNt 1.0, 2, 0,0! 

CALL  CST«INGI280,973, 37hCCRRELATlCN  FUNCTION  OF  THE  PROCESSi.l 
CALL  CSTRInG(496,7S,  ICHT IKE  LAGS.! 

CALL  MAPS(0.,25.,-1.5,1.S,.15,.^5,.IS,.9I 
CALL  LTNEIO.,0.,23.,0.) 

XOLO-0. 

yoLO>o. 

DO  I 1^1, INOOUT 

1 SUM»SUM*OUTPUT(  I 1 
KN>  INOOUT 
AVE*SUM/RN 

00  2 1=1,26 

IMl»I-l 

SUH=0. 

ISTCP»IN00UT-IH1 

00  3 rT=i,isrop 

IN0X,,IT*IM1 

3 SUM=SUM*( OUTPUT  « I T 1~AVE !•( OUTPUT ( INCX !-AVE! 

SPECT»SUM/RN 

IFIIKl.EC.O!  SP2ER0*SPECT 

SPECT=SPECT/SP2£R0 

RX-IMl 

CALL  LYKE(XCLO,YCLO,RX,SPECT) 

XOLO=RX 

VOLCxSPECT 

2 kRlTE{6,5!  IMl.SPECT 

5 FORMAT  (I5,F20.B! 

WRITEf6,7!  SPZERC 

7 FORMATKO  THE  VARIANCE  OF  THIS  SAMPLE  HAS  *,F15.8I 
CALL  OPTION  (1,0,2,1,01 

CALL  CSTRING(75,360,25HCCRRELAriON  COEFFICIENTS.) 

CALL  FRAME 

RETURN 

END 


1 


